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ABSTRACT 


The initial conditions (IC) used in Numerical Weather 
Prediction (NVP) codes have a degree of uncertainty, due to 
observational and other errors. Given this, it is futile and 
computationally expensive to increase the resolution of these 
numerical simulations beyond the point where the resolution errors 
becomes secondary to initial condition dependent (ICD) errors. In 
an attempt to find an optimum resolution, in terms of a spectral 
truncation number N, the resolution and ICD errors in a Spec^iral 
Barotropic Vorticity model is studied. A case is made as to why 
the study of this simple model yields the upper limit of the 
resolution required by more realistic models. It is found that, 
given error in the IC of 5-10% RMSCroot mean square), N=32 
produces more than adequate resolution for three day predictions 
of the velocity field. The vorticity field, however, requires 
higher resolution. It is found that the the errors induced by 
truncation of IC grows quite slowly with time. Another interesting 
observation is that, for the same overall RMS error, the 
predictions are much more sensitive to errors in the lower order 
spectral coefficients than in the higher order ones in the initial 
conditions . 



CHAPTER ONE 


INTRODUCTION 


Historical Background : 

The cornerstone of the scientific weather forecast is tl 
so-called synaptic (same-time) map of the current conditions (i.e 
pressure, velocity etc.) over an area spanning distances of tli 
order of a thousand kilometers (called the synoptic sca.ls:> in eac 
direction. These maps can then be used to predict the weather fo 
the following few days, by the use of either the forecaster’s skil 
which combines experience, scientific reasoning and 
we 1 1 -developed intuition, or by the use of numerical models. 

However, the first synoptic map of surface pressure was no 
constructed until 1820 and then only with data nearly 50 years old 
It was only after the invention of the Telegraph in 1845 tha 
synoptic maps could be constructed fast enough to be used i 
operational forecasting. During the period called the «mpCirica.L er 
(1850-1920), weather organizations to provide forecasts for th 
general public, agriculture, flood warning, transportation, etc 
were founded. But the forecasting in this period relied solely oi 
the skill of forecasters in interpreting synoptic charts and othei 
observed data. 

Bjerknes, in 1904, realized that the motion of the atmospheri 
is fundamentally a solution of the basic equations of fluic 
mechanics. Bjerknes also realized that this system of non-lineai 
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partial differential equations did not posses a general solution, 
and that, in addition, the data-gather ing capability in his time 
was wholly inadequate to determine the initial conditions for the 
initial value problem involving the fluid flow equations. 

Less than two decades later, L. F. Richardson attempted to 
solve the system of equations numerically using manual calculations 
(which involved months of tedious work). This pioneering effort, 
published in 1922, unfortunately predicted pressure changes of an 
order of magnitude greater than actually observed. It is now known 
that Richardson had made a fundamental error in his approach - the 
numerical time-stepping scheme used was ii-nstaJbi& (fundamental work 
on numerical instability was done by Courant, Friedrichs and Lewy a 
few years after Richardson had published his work). But the result 
of the first numerical simulation of the atmosphere was to 
discourage other similar attempts for three decades till the 
invention of the first practical electronic computers. However, the 
invention of the radiosonde, during this period, led to further 
empirical improvements in forecasting. Meanwhile, theoretical 
research laid the foundation for a departure from pure empiricism. 

Electronic computers provided the last ingredient needed for a 
breakthrough. In 1950, Charney, Fj0rtoff, and Von Neumann produced 
the first successful numerical forecast using an 
equival ent-barotropic vorticity model which was far simpler than 
the one developed bv Richardson. The weather predictions made with 
this simple model, however, were nearly as good as those made by 
highly skilled forecasters using empirical techniques. With this 
success, numerical weather prediction (N¥P) was soon on the way to 
becoming the primary tool for modern weather forecasting. 
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During the last four decades there have been rapid advances in 
all aspects of NWP , supported by almost incredible developments in 
the speed and memory size of electronic computers. The development 
of meteorological satellites has provided global observations that 
have spurred new insights about the weather. Satellites have also 
been very valuable for short-range disaster warnings of hurricanes 
etc., especially over and near areas with few other observational 
platforms. However, data-gather ing from these satellites has not 
yet reached the sophistication and accuracy needed in providing 
initial conditions for NWP. 

Numerical weather forecasting requires a closed set of 
governing equations (which are generally partial differential 
equations), suitable initial and boundary conditions, and a 
numerical method for solving the initial boundary value problem 
implied in forecasting. The boundary conditions provided in a 
global model will be the natural periodic boundary conditions. 
However, if NWP is being attempted using a limited area model, some 
artificial boundary conditions are devised so as to be able to 
solve the problem. 

Obtaining the initial conditions for the numerical prediction 
codes is not a trivial task. It involves the collection and 
checking of surface- and upper-air meteorological observations, 
which work forms the bread-and-butter activity of most 
meteorological stations. However, this data cannot be directly used 
as initial conditions because, firstly, the initial conditions are 
to be specified over a given set of grid points covering the 
forecast region, which generally will not coincide with the 
positions of the observation stations. Therefore, some form of 
interpolation has to be used. Furthermore, the raw or interpolated 
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data, due to observational or interpolation errors, may contain 
energy in certain scales (length or time) which cause the initial 
data to be unphysical or which may cause the computations to become 
unstable. Therefore the raw data is subjected to dyuaraical and 
statistical constraints, which are used to ’massage’ the data into 
a form in which the unphysicai or dangerous modes are removed. 


Sensitivity to initial conditions: 

A fundamental and far-reaching property of the equations 
governing atmospheric motion is the sensitivity of their solutions 
to small perturbations in the initial state. Solutions initially 
very close to each other evolve so that the difference between them 
is as large as the difference between randomly chosen states of the 
system . 

There is motion in the atmosphere on scales ranging from 
planetary waves of thousands of kilometers to mi 1 1 imeter-soale 
turbulent eddies. The initial state for the numerical weather 
prediction models, however, is obtained from observations of only 
those scales of motion large enough to be resolved by the observing 
network. Smaller scales of motion are completely ignored and 
indeed, through alias error. cause uncertainty in the 
determination of large scales . Therefore, no matter how accurately 
a numerical model computes the evolution of the atmospheric motion 
there remains a fundamental uncertainty in the forecast due to the 
error in the initial state. 
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Predictability of the atmosphere: 


It has long been realised that the equations of the atmosphere 
are highly sensitive to perturbations. Small differences in the 
initial state tend to magnify and finally lead to very dissimilar 
final states after a period of a few days or at most a couple of 
weeks. This has lead researchers to speculate about and seek the 
limits of predictability of the atmosphere. This research is aimed, 
in part, at finding a time-scale over which even the smallest 
errors in the initial state would lead to significantly different 
final states. 

Thoiripso-n [1] pioneered the study of predictability by 
examining the behavior of an idealized two dimensional barotropic 
atmosphere for an initial error that is isotropic and homogeneous. 
Loj'tf-nat [21, using a turbulence closure scheme for two-dimensional 
flow, showed that error initially confined to smallest scales of 
motion (order of a few meters) ’cascaded’ to higher and higher 
scales, inducing significant errors in the so called "cumulus 
scales" (2—10 km) within half an hour, in the "synoptic scales" in 
two days, and in the planetary scales after two weeks. 

Kraichnan [31, Leith [4] and Leith and Kraichnan [5] did 
similar calculations using sophisticated closure schemes. These 
studies concluded that, in three dimensions, whatever the scale of 
initial error, the error will grow and finally dominate the large 
scales of the flow. Thus a three dimensional flow field is not 
"predictable" for all time even if the initial field is known to an 
arbitrary degree of accuracy. 
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Studies which use models closer to the real atmosphere have 
also been done to discern the effect of chosen atmospheric factors 
on predictability. The effect of Rossby waves, and the Coriolis 
parameter ('beta'), on error growth in barotropic models was 
studied by Sasctevant et al , [61 and by Hallcjway [7]. Both studies 
agreed that Rossby wave propagation enhances predictability (, i . e . 
slows error growth). Vallis 18] found that such increased 
predictability due to Rossby waves is not found in baroclinic 
mode 1 s . 

IZoaxis [9] used a two-level quasi-geostrophic model to study 
predictability in the extended range. He found that it took about 
30 days before the forecast skill of the model was about the same 
as a persistence forecast. Strcruss flOl also used a two-level 
quasi-geostropic model to study the effect of baroclinic 
instability and wave-wave interactions on the predictability. 

Similar predictability studies have also been attempted using 
more complex General Circulation Model s(GCM). Charne?y t?t al . [11], 
SmagoTiTvshy [12] and Williamson and KasaHara, [13], attempted GCM 
runs with initial conditions that were in some sense close to each 
other. The doubling time for small errors was found to be quite 
model dependent. The role of moisture was included in the later GCM 
studies of Ba^ds »t a.1 . [14]. 


Errors in the initial condition : 


The studies mentioned above 
errors in forecasting i.e. the 
magnify small differences in the 


concentrate on 

one 

cause of 

the 

tendency of 

the 

atmosphere 

to 

initial state. 

However, even in 

a 
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perfect model, the accuracy of a numerical forecast will depend on 
at least two factors, namely, (a) the growth rate of error as 
determined by the equations of motion, and (b) the magnitude and 
scale of the initial error as determined by the accuracy and 
resolution of the observing network. The first of these sources of 
error have been studied by those interested in the ‘predictability’ 
of the atmosphere. ¥e shall now consider the second. 

The process of converting observed data to the initial 
conditions for a NWP code is essentially one of finding a suitable 
analytic representation for point-wise data. It is a process of 
curve-fitting subject to dynamical and statistical constraints. 
Often, a least-squares type approach is used to achieve this end. 
However, there are many sources of error on the path from the 
‘real* initial fields to the initial conditions finally used in the 
numer i ca 1 mode 1 . 

To get a good interpolation function for point-wise data, it 
is usually necessary to have a systematic distribution of the 
data-points. However meteorological observation stations have a 
very uneven distribution over the earth, with large land masses in 
the world, and the oceans which cover three-fourths of the globe, 
scarcely represented. Interpolation using unevenly scattered 
data-points leads to errors in scale, by which energy at a certain 
length scale is mistakenly attributed to another scale. 

Because the weather stations are at least of the order of 
hundreds of kilometers apart, only the large scale structure of the 
atmosphere may even be attempted to be captured by observations. 
However, large-scale motions may only be detected accurately by 
observations over correspondingly large time-scales. This however 
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is not generally feasible, especially in the upper-air 
observations. Therefore, small time-scale behavior may be mistaken 
for large time-scale behavior, thereby compounding aliasing error 
problems. The small-scale, high-wave-number , structure of the 
atmosphere cannot be obtained due to the poor resolution of the 
observation net-work; and hence a not insignificant portion of the 
kinetic energy may be ignored even in the best case. 

Observation data undergo a detailed processing, which involves 
checking the location of observations, for gross errors, for 
hydrostatic consistency, etc. After preliminary error checks, the 
fields are subjected to ‘ Objectiz^e- Analysis'. The aim of objective 
analysis is to extract the maximum information from the 
observations, consistent with climatological records, space 
correlations between the meteorological variables, etc. This 
requires a knowledge of the statistical structure of the fields of 
the variables, which may be analyzed separately as in -unixiarrate 
analysis, or together as in miul ttvariatef objective analysis. 


However, objective analysis procedures are generally not 
enough to provide the initial conditions for an N¥P code. In the 
atmosphere, the wind and pressure fields are in a 
quasi-geostrophic, quas i-nondivergent state and little energy is 
associated with the ageostrophic motion. However, the objective 
analysis schemes generally do not yield mass and motion fields that 
are in the delicate state of balance that exists in the 
atmosphere. Consequently, the use of only objectively analyzed data 
to initialize an NWP model may generate large, spurious 
inertial-gravity modes, if the model permits such modes. 
Furthermore, because of numerical errors, computational 
instabilities , etc , there may be spurious growth of short waves, 
which if left unchecked could ’blow up’. Even otherwise. it is 
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usually desirable to eliminate very high wave number components. 
Therefore, the fields have to undergo a further process of 
’initialization’ which will eliminate these unphysicai and 
troublesome aspects of the objectively analyzed data. 

It may be remarked that the above procedure for generating 
initial conditions does not actually correct errors in the raw data 
(other than in locating gross errors). What it actually does is 
generate fields which are similar to the observed ones but which 
have dynamic and statistical consistency and which will not cause 
numerical difficulties. Therefore it may be argued that the 
initialized fields are far from being error-free. 


Initial conditions errors versus numerical integration errors: 

Since the 1970 ’s there has been a shift in the numerical 
method used for global atmospheric models from the 
finite-difference to the spectral method. The spectral method has a 
number of advantages, primary amongst which is accuracy. Th3s 
method has the so-called ’exponential’ or ‘spectral’ accuracy which 
ensures that, for smooth solutions, the error in the spectral 
representation decreases faster than any finite power of the 
grid-interval. Thus, as the number of grid-points are increased in 
the numerical scheme, the accuracy of the representation rises 
dramatically. Furthermore, in global spectral models, generally a 
Spherical Harmonics representation is used which ensures that there 
are no problems of the discontinuity of the dependent variables and 
their derivatives at the poles. 

The principal disadvantage of the spectral representation is 
that a large amount of arithmetic is necessary to calculate 
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spectral equations with even a moderate number of retained modes. 
The transform method described by Or^siaj^l i 5 ] reduces the number of 
arithmetic operations (per level, per time-step) from O(N^) to 
0(N*) where N is the spectral truncation number. However, as these 
many operations have to be done every time-step, increasing N, 
which undoubtedly increases accuracy rapidly, also extracts a 
progressively heavy price in the number of computations. 

With the developement of these highly accurate numerical 
methods, a certain imbalance has crept into the use of numerical 
weather prediction codes. It is much easier, given the rapid rate 
of increase of the speed and power of the computers, to increase 
the formal accuracy of the numerical integrations performed by the 
codes, by increasing the number of grid-points, than to increase 
the accuracy of the initial conditions being fed into them. This 
increase in the formal accuracy, of course, increases the 
computational burden. The result is that the computational 
requirements of NWP codes grow to fill the capacity of the latest 
and most powerful computers. We may ask : to what purpose? Does 
improved accuracy of the numerical integration necessarily increase 
the accuracy of the forecast, given the uncertainty in the initial 
conditions ? This study seeks a partial answer to this question. 

The error in a forecast may be divided into four categories 
which are mutually independent, i.e, (a) errors due to 
faults in the model, say, by incorrect representation of the 
physics, (b) timt? errors due to inaccuracies in the 
temporal integration of the model equations, (c) r^ssol^it ion orx'orsi 
caused by inaccurate representation ol the spatial fields, and 
finally (d) initial condition d&pe>ndf=nt (ICD) errors caused by 
incorrect specification of the initial conditions to the problem. 
We are primarily concerned with the last two types of error in this 
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study. The aim of this work is to investigate tlie degree of 
spectral accuracy that would be sufficient to resolve a numerical 
simulation of the atmosphere with a specified amount of error in 
the initial conditions. 

We are primarily concerned with a comparison of the resolution 
errors with the initial condition dependent errors. Say, with a 
given number of spectral coefficients CN^, say) we. obtain 
resolution errors which are approximately the same as the expected 
ICD errors then, because of the exponential accuracy of the 
spectral method, we can decrease the resolution errors 
substantially by a moderate increase in N. But as this will only 
increase resolution accuracy, while leaving the initial condition 
dependent (and other) error essentially unaffected, it would be 
pointless to increase N too much; this will only increase the 
computations for the illusion of greater accuracy. Therefore the 
optimum value of N, which we will call N , will be slightly 

Opt 

larger than that value of N at which resolution and ICD errors will 
be comparable. 

As this is a numerical study, we need a numerical model that 
suits our needs. We have chosen to study simulations of the 
one-level barotropic vorticity equation. There are two reasons for 
this choice: 

Firstly, the model is a very simple one and requires the least 
amount of computation for a given number of spectral coefficients. 
The second reason, which is more compelling, involves the intimate 
relationship between the predictability of the atmosphere and 
initial condition dependent (ICD) error. It is clear that for a 
less predictable model atmosphere there will be a greater rate of 
divergence of two states starting from almost identical initial 



conditions. This means that the same error in the initial 
conditions will cause greater TCD error in a less predictable model 
atmosphere than in a more predi cable one and ui ce-ner-va. 

Now, the studies of predictability quoted above have 
established that (a) three dimensional atmospheres are less 
predictable than two-dimensional ones, and (b) introduction of 
barocl inicity makes the model atmosphere less predictable. ¥e may 
conclude that the one-level barotropic model is mor© predictable 
than the real atmosphere (and more realistic models), and will, 
therefore, have a less rapid growth of ICD errors. However, the 
resolution errors in this model will be about the same (for the 
same number of spectral coefficients per level) as in any other 
model. In comparing the resolution and ICD errors in a barotropic 
model we will, therefore, be ■undLer^stiTna.tiTis the relative magnitude 
of the latter (with reference to a more realistic model of the 
atmosphere). Therefore the optimum resolution N ^ , which is 
chosen such that the resolution errors are dominated by the ICD, 
will be an ovsi'tgsiimai^. Therefore, in using the barotropic model 
we are finding the -upper limit of the spectral resolution required 
in N¥P codes. 

The approach that is followed in this study is similar to the 
’identical twin’ numerical exoeriments that are so common in 
predictability studies. First, certain ’standard’ initial 
conditions are selected. The highest resolution calculations of the 
evolution of the model atmosphere, using these standard initial 
conditions, are treated as the ’exact solution*. Other simulations, 
done using either (a) coarser grids (poorer resolution) or (b) 
erroneous initial conditions, are compared with the exact solution 
to determine resolution errors and ICD errors, respectively. Error 
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in the initial conditions is introduced by randomly perturbing the 
initial spectral coefficients from the ’standard’ values. Uniformly 
distributed random numbers are used to generate errors of any 
chosen standard deviation. In other simulations, the error is 
introduced in low, and high, wavenumber initial spectral 

coefficients, respectively, to assess the wavenumber dependence of 
the rate of growth of the initial error. 

The barotropic model is one of the oldest numerical models of 
the atmosphere and has been studied extensively. Detailed 
explanation of the spectral form of this equation and the numerical 
technique used in solving it are given, for example, in Hal i in&z- 
and Williams [16], Holton [17], and in Washington and 
PoT-hinsonl 18] . 

The thesis contains a detailed discussion of the derivation of 
the nondivergent barotropic vorticity equation, an introduction to 
the spectral method, spherical harmonics, and the energy spectra 
of scalar fields in chapter 2. Chapter 3 contains a description of 
the non-dimens ional ization of the barotropic vorticity equation, 
its spectral formulation , and the algorithm used in the present 
work. Finally, chapter 4 lists out the results and related 
discussions . 
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CHAPTER TWO 


THE METHOD 


This chapter contains the details of the methods implemented 
in the present work. The topics covered include the derivation of 
the synoptic scale momentum equation in spherical coordinates, the 
Barotropic vorticity equation, the discrete Fourier transform. 
Associated Legendre polynomials, Gauss lan-quadrature , Aliasing, and 
Spherical Harmonics. 

C2.1> Primitive Equations 

Atmospheric motions are governed by the principles of 
conservation of mass, momentum, and energy. The momentum equation 
relates the rate of change of momentum (i.e. the total time 
derivative of momentum) in an inertial reference frame to the sum 
of the forces acting on the fluid. By the total time derivative of 
a quantity (also called the Lagrangian derivative) is meant the 
rate of change of that quantity as seen by a particle of the 
moving flvicL. This is related to the time and space derivatives as 
seen by a stationary observer in the following manner; 

d a ^ ^ ^ ^ ^ ^ 
dt at ax a y a z 

where u, v, and w are the component velocities in x, y, and z 
directions respectively. The quantity on the left side of the 
equation is the total derivative; the first quantity on the right 
is called the local derivative and the last three terms are 
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is called the local derivative and the last three terms are 
called the convection or advection terms (which are often denot-ed 

A 

by the notation V- V ). 


In meteorology it is desirable to describe the motion in a 
reference frame rotating with the earth. The transformation of the 

A 

total time derivative of a vector A in an inertial frame and the 
corresponding total derivative in a rotating system with constant 
angular velocity n is 


dA 1 

d t J tri&rtval 





d t J rotating 


A. A 

+ O X A 


( 2 . 1 . 1 ) 


In an inertial frame Newtons second law of motion may be 
written aymbolically as 


( 2 . 1 . 2 ) 

The left hand side represents the rate of change of the absolute 

A 

velocity U following the motion as viewed in an inertial 
system. The right-hand side represents the sum of the real forces 

A 

acting per unit mass. By applying (2.1.1) to the position vector r 
for an air parcel on the rotating earth we get : 


If 1 . + n X r (2.1.3) 

dt Jtnamal dt J rotating 

The first term on the RHS of (2.1.3) is the velocity of the 
particle with respect to the rotating system which will be denoted 

A 

by Ua. We can then rewrite (2.1.3) as 


Ua = U + n X r (2.1.4) 
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which states simply that the absolute velocity of an object on the 
rotating earth is equal to its velocity relative to the earth plus 
the velocity due to rotation of the earth. 


Next we apply C2.1.1) to the velocity vector Ua and obtain 


-j - i— . . , = 5-T- I + Q X Ua 

at JvtxarUal at J rotating 


Substituting from (2.1.4) into the right-hand side of (2.1.5) 
we get 


dUo 


5 


dT j vr^rtLal 


-dCU ♦ O X r) 
d t 


+ n X (Uct X r ) 


^ + 20xU-0*R 


( 2 . 1 . 6 ) 


where the last expression on the S.H.S. is in the rotating system. 
Here R is the projection of r on a plane perpendicular to the axis 
of rotation, Therefore the acceleration of a body with respect to 
an inertial system equals its acceleration in a rotating system 
plus the Coriolis and centripetal acceleration terms. The second 
and third term on the RHS of (2.1.6) are called coriolis and 
centripetal acceleration respectively. These are called 
’fictitious’ acceleration terms because their existence is solely 
due to the fact that motion is being observed from a non-inert lal 
frame . 


If we assume that the only ’real’ forces acting on the 
atmosphere are the pressure gradient force, gravitation, and 
friction, we can rewrite Newtons second law (2.1.2) as 


IG 



(2.1 .7) 


dU 


2 n X u 


(1 / p) V P + g + Fr 


where Fr designates the friction force, P is the pressure and the 
centripetal force has been combined with gravitation in the gravity 
term g. Equation (2.1.7) is the fundamental form of momentum 
equation dynamic meteorology. 


(2.2) The Component equations in Spherical coordinates 

The shape of the earth is nearly spherical, therefore it is 
convenient to expand (2.1.7) in spherical coordinates. The 
coordinates are then CX,0,Z), where x is the longitude, 4 > is the 
latitude and Z is the vertical distance above the surface of the 
earth. If the position of amibject is denoted by (\,^,Z) and if 

^ Jk. Jk. 

the unit vectors i, j ,k are now taken to be directed eastward 
(increasing X), northward (increasing ^) , and upward (increasing 
Z) , respectively, the relative velocity of the object with respect 
to the earth is 

^ /K A. Ay 

U»iu + jv + kw 

where the components u, v, and w are defined as 

u s r oosd> , V = r , w s -^| — (2.2.1) 

Bere , r is the distance to the center of the earth which is 
related to z by r = a+ z, where a is radius of the earth, 

jenerally, the variable r in (2.2.1) is replaced by the constant 

JL. This is a very good approximation since z « a for the regions 
if the atmosphere with which meteorologis t is concerned . For 
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convenience, meteorologists have traditionally uffed a coordinate 


system which is spherical in nature but which is cartesian-l ike 
over small distances. This is done by defining a coordinate x as 
the eastward distance and y as the northward distance from a given 
point (dx= a cos<^ dX, dy= a d^) . For small areas, as the curvature 
of the earth has negligible effect, the coordinates x and y 
essentially define a locally cartesian system. The velocity 
components can be expressed as u s (dx/ dt), v = C dy/ dt), and w 

A A 

s ( dz/ dt ) and the unit vectors i, j, k are in the x, y, and z 
directions, respectively. However the system is not cartesian as 

A. A A 

i, j, and k depend on the location of the point (x, y, z). This 
position dependence of the unit vectors roust be taken into 
account when the acceleration vector is expanded into its 
components on the sphere. Thus we write 


.A* 

dU 

“arr 


= i 


du 

dt 


+ j 


dv 

dt 


+ k 


dw 

dt 


u 


An. 

di 

dt 


dj 

dt 


+ w 


As 

dk 


dt 

( 2 . 2 . 2 ) 


where 


di 

dt 

A 


= U 


d i 
9 X 


= u 


a cos^ 

A 


ii- = u 

5 j 

-> Hh 

V ^ J 

dt^ 

9 X 


5 y 

A 

dk _ 

u 

A 

i 

+ JL. 

St 

a 

a 


(j sin^ - k cos^) 


Substituting equation 
components on the left 
equation we get 


( 2.2.2) into 

and right hand 


(2.1.7) 
side of 


and 

the 


(2.2.2a) 

(2.2.2b) 

(2.2.2c) 

eK}uat ing 
resulting 


du 

Wl 


u V tan<;i> ^ u w 
a a 


£p 
p dx 


2 Q w GOS0 +20 V sin^ +Fx 

(2.2.3) 
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dv 

dt 


2 Cl u 8in<j6 +Fy 


(2.2.4) 


. u* tan^ . vw 

T T 

a a 


1 £E 

p ^ 


dw _ 
dt a 


= — ^ - g + 2 O u ooBtp +F 2 
p&z 


(2.2.5) 


where Fx, Fy, and F* are the three components of the frictional 
force. The terms proportional to 1/a on the LHS in (2 . 2 . 3)-( 2 . 2 . 5 ) 
(which are called the curvature terms because they arise due to 
curvature of the earth) are non-linear. Although, as will be shown 
in the next section, the curvature terms are unimportant for mid 
latitude synoptic scale motions and may be neglected, 
(2 . 2 . 3)-(2 . 2 . 5) are still non-linear PDEs as can be seen by 
expanding the total derivatives into their local and advective 
parts. In general the non-linear advective acceleration terms are 
comparable in magnitude to the linear local acceleration terms. 
It is primarily the presence of non-linear advection processes 
that causes the rich and complex behavior of the atmosphere. 


(2.3) Scale analysis of Equation of Motion 

The equations of motion derived in the last section are 
complicated coupled non-linear PDE’s. A natural hope that arises is 
that a way may be found to simplify them. This may be done by 
evaluating the order of magnitude of the various terms in the 
equations and neglecting those that are small in comparison to 
others. This approach has been found to be very fruitful in 
meteorology and often leads to a substantial reduction in the 
number of terms in the equations. This kind of order of magnitude 
inalysis is called scale analvsis or scaliTig. 
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In scaling, typical values of the magnitude of the field 
variables, and the amplitude of their fluctuations, and the 
characteristic length, depth and time-scales are estimated. These 
typical values are then used to determine the magnitude of various 
terms in the governing equations. 

The character of the atmospheric motions depends very 
strongly on the horizontal scale. Table (2.1) presents various 
types of motions classified by horizontal scale for the spectral 
region from 10"’ to lO^m. 


TABLE- 2.1 


Types of motion 

Horizontal scaleCm) 

Kolecular mean free path 

10"'' 

Minute turbulent eddies 

10"*-10"* 

Small eddies 

lO' ^-1 

Dust devils 

1-10 

Gusts 

10-10* 

Tornadoes 

10* 

Cumulonimbus 

10® 

Fronts 

io'‘-io'’ 

Hurricanes 

10® 

Synoptic cyclones 

10® 

Planetary waves 

10* 


Elimination of terms by scaling not only simplify the 
lathematics, but in some cases it completely eliminates, dt* 
ilters, unwanted type of motions. For example, in an analysis of 
urricanes (~10° m) we are usually not interested in motion of the 
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)rder of lO^m (gusts) and these may be filtered out of the 
equations by scaling. 

The term synoptic is associated with the analysis of 
>bservations taken over a wide area at or near the same time. The 
lynoptic scale is the characteristic scale of disturbances which 
,re depicted on weather maps(10**m). In order to simplify 
2.2.3)-(2.2.5) for synoptic scale motions we define the following 
ypical values of variables based on observed values for 
id latitude synoptic systems. 

U ~ 10 m e”* horizontal velocity scale 
W ~ 1 CTO e”* vertical velocity scale 
L ~ 10^ TO length scale [ ~ 1/ (2 rr) wavelength 1 

D ^ 10'* w depth scale 

Ap/p ~ 10® TO® e”® horizontal pressure fluctuation scale 
L/U ~ 10** e time scale 

j take a typical value of the latitude as 45° and estimate 

le Coriolis parameter fs 2Q sin^ at lO”* s”*. Now we can estimate 
le magnitude of each term in the equation (2.2.3)-C2.2.5) by 
‘.placing each variable by its typical value and replacing 
>rivatives by the ratio of the typical values of the 
fferentiated variable and the differentiating variable. Table 

1.2) shows the typical magnitude of each term in (2.2.3) and 

1.2.5) obtained by this procedure. 

Frictional terms are neglected because, for synoptic time 
ales, frictional dissipation is significant only in the region 
ry near the surface of the earth. From table(2.2) it is clear 
at the coriolis force and the pressure gradient force are of the 
eatest magnitude while the other terms are insignificant in 
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comparison. Therefore, retaining only these two terms 
and (2.2.4), we obtain as a first approximation the 
re lat lonship 


in (2.2.3) 
g&os i roph. t c 


f V 


P 


£p 

dx 


f U 2* 


1 £P 
P ^ 


(2.3.1) 


where f = 2Do sin^ is called Coriolis parameter. The geostrophic 
expression gives the approximate relationship between the pressure 
field and horizontal velocity in synoptic scale systems. However 
it contains no reference to time, and therefore cannot be used to 
predict the evolution of the velocity field. It for this reason 
that geostrophic relationship is called a diagnostic relationship. 


To obtain a prediction equation it is necessary to retain the 
acceleration term in (2.2.8) and (2.2.9) .The resulting 
approximate horizontal momentum equations are 


TABLE- 2.2 Scale-analysis of Horizontal Momentum Equations 


X and y momentum equations are : 


du 

u V tan^ 1 
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1 ^ - o 

fie W COSip ^ 
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Ob V 

sin(^ 

dt 
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a 

p ex 
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+ 

tan^ 
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_ 1 dp 
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Qo u 

sin^ 

dt 

Scale of 

a 

individual 

a 

terms : 

" p dy 



U*/ L 

U*/a 

UV/a 

AP/(pL) 





Magnitude of terms ( 

ms ) : 






10'* 

10'® 

lO'* 

10"® 

10"® 


10" 

•3 
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du 


f V 


1 ^ 
p dx 


(2.3.2) 


dv 

dt 


+ 


f u 


1 £e 

p ^ 


(2.3.3) 


As table (2.2) shows, the acceleration terms are about an order of 
■magnitude smaller than the coriolis force and pressure gradient 
force. Thus, a small error in measurement of either velocity or 
pressure gradient will lead to very large error in estimating the 
acceleration by using the above equations. 


(2.4) Barotropic Vorticity Equation 

The equations (2.3.2) and (2.3.3) may be used to develop a 
very simple model of the atmosphere called the Barofropic model. 
This model assumes that the velocity field above a point (x,y) can 
be represented by a single value, i.e. it ignores the vertical 
variation of the vertical variation of the velocity field and 
approximates a 3-D field by a 2-D one. A further assumption in the 
model is that the 2-D velocity field is incompressible and 
frictionless. Physically, these assumptions equate the atmosphere 
to an inviscid barotropic fluid i.e. one in which density is a 
function of pressure alone (and independent of say, temperature). 
It has been found that at any given time, significant portions of 
the atmosphere indeed exhibit nearly barotropic behavior. The model 
is, therefore, not as unphysical as it may at first appear. 

Flow fields can be decomposed into the sum of two vectors: one 
nondivergent( incompressible) but with vorticity, the other 
divergent(compressible) but irrotational . For large-scale weather 
patterns in mid-latitudes, the nondivergent , rotational vector is 
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usually by far the larger of the two, so that by studying vorticity 
we are studying a representation of most of the flow. It is well 
known that the motion of a 2-D incompressible fluid can be 
represented by a single scalar function called the stream function. 
The stream function, by its very existence, ensures that the 
velocity field described, complies with the incompressibility 
condition. Therefore, there is no need to solve the continuity 
equation if a stream function formulation is used. The two n»omenttHB 
equations for u and v velocity are replaced by two scalar equations 
relating stream function and vorticity which will provide an 
equivalent description of the motion. ( Vorticity is a vector for a 
3-D field but in 2-D it has only one non-zero component and can 
hence be regarded as a scalar ). The velocity u and v can be easily 
obtained from the stream function if needed. 


The vorticity equation can be derived using the momentum 
equations (2.3.2) and (2.3.3) . We differentiate the x component 
equation with respect to y and y component equation with respect 
to x : 


g f g u 

dy ^ d t 


u 


+ V 


d u 

Q y 


+ W 


d u 
9 z 


f V 


1 £_p \ 

P & X J 


(2.4.1) 


9 f a V 

'Sr\~S~T 
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a V 

■rr 
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a V 

a y 


w 


a V 
a z 


f V 


1 t-p \ 
p ^ y / 


(2.4.2) 


Subtracting (2.4.1) from (2.4.2) 
component of) vorticity C = ( ^ v/ 
the vorticity equation [ using 
parameter depends only on y so that 


and defining the (vertical 
dx)-(^u/dy)we obtain 
the fact that the Coriolis 
df/dt = V ( df/dy ) ] 
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d(C + f) 


) 


( 


) 


= - (C + f) ( 


d u , d V 

TTiF Try 


& w d V _ 0 w 0 u 
0 X d 'z ~ ^ y a z 


+ (1/p^) ( 


a p a p 

W~x T-y 


a p a p 

3-y anr 


) 


(2.4.3. 


The three terras on the right of (2.4.3) are called the 
term, the tilting term or the twisting term , and the 
term, respectively. 


divergence 

solenoidal 


Equation (2.4.3) can also be simplified using scale analysis. 
For synoptic scale motions we can choose the same characteristic 
scales as given in section (2.3). 

U ~ 10 we”* horizontal velocity scale 

W ~ 1 ciY. s“* vertical velocity scale 

L ~ 10** w length scale [ ~ 1/ (2 n) wavelength] 

D ~ 10* w depth scale 

ap ~ 1 kPa horizontal pressure fluctuations 

p ~ 1 kg mean density 

ap/p ~ 10"* fractional density fluctuation 

Ap/p ~ 10* Tn*o”* horizontal pressure fluctuation scale 

3 

L/U ~ 10 a time scale 

f^ ~ 10 Coriolis parameter 

df/dy * /3 ~ 10~** "beta" parameter 


Using these scales to evaluate the magnitude of the terras in 
(2.4.3) we note that 


C ~ 10"° s"‘ 


f ( 


0 U ^ 0 V 
ax ay 


- 10"*’ 
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10 


-11 


Retaining now only the term of order 10“*° s“* in the vorticity 
equation, we obtain as the first approximation for synoptic scale 
mot ions 


dCC + f) _ 

dTT 


(C + f) 


( 


3 u d V 

3 X 3 y 


(2.4.4) 


As a first approximation the change of absolute vorticity 
following the horizontal motion on the synoptic scale is entirely 
due to divergence effect . 


BAROTROPIC MODEL : 


We use the one parameter Barotropic model for our study. By 
one parameter model is meant that the model uses only one data 
point in the vertical direction. The assump^tions in the Bar-otropic 
modal are 

(i) Inviscid and incompressible flow. 

(*) Motion is two-dimensional (vertical motion is 
altogether neglected). 

(a) Effect of frictional forces is neglected. 

Then vorticity equation becomes 

d(C ■«- f) _ n 
d t " ^ 

or I = - U . V ( C + f ) (2.4.5) 
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Note U in the case of barotropio flow will consists of only 
nondivergent part of the velocity field. Furthermore, as only 
one level is allowed , we may represent the velocity field at 
some chosen isobarCsay 300 or 500 mb) where the field is nearly 
divergence less . 

The barotropic model is a useful approximate forecast equation 
for short range weather predictions. To a limited extent model is 
used to understand the basic idealized properties of large scale 
dynamical motions. 

(2.5) Spectral Method 

Spectral method are global approximations unlike, finite 
difference methods which make approximations of, say, the derivative 
of a function by using the nearest neighbors. The spectral approach 
is to fit a function through all the grid points and use this to 
find the derivative at any point. 

The essence of the spectral method lies in the approximation 
of functions by a truncated series of standard known functions 
called trial functions. The expansion of a function A(x) would be 
of the form 

A Cx) = E A Z (X) , a S X b (2.5.1) 
app N 

where A (x) is the spectral approximation to A(x) and A^’a are 
(complex) constants called spectral coefficients. Z^(x) is a trial 
function. The wave number n varies from 0 to N-1 tor from (N/2)-l 
to N/2 depending on the type of trial function] where N is the 
total number of terms in the series. 
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The commonly used trial functions are 

(a) Fourier 

(b) Chebyshev 
Cc) Legendre 

(d) Spherical harmonics 
Ce) Hermite 
(f) Laguerre . 

The trial functions are orthogonal functions. Not all the 

trial function sets will give infinite order accuracy but the ones 

named above will do so. By "Infinite order accuracy" is meant that 

the error in the expansion will fall faster than any any 

polynomial of (1/N) as N — >cd. It is also called "exponential" and 

"Spectral" accuracy. The choice of the sets of trial functions 

depends on the problem. Fourier functions can only be used for 

expansion of periodic functions, while spherical harmonics can 

only be used for the expansions in a spherical coordinate system. 

Chebyshev functions are suitable for finite domain (rectangular) 

problems. The choice of function set is also dependent on the 

availability of efficient algorithm to evaluate spectral 

coefficients A for any given function A(x). 
n 


(2. 5. a) Evaluation of coefficients 

Obviously, the crux of finding good approximations for A(x) 

in equation (2.5.1) lies in finding the the best values of ( the 

as-yet-unknown) spectral coefficients. This is where the 

orthogonal property of the trial functions are useful. ¥e can 

utilize the orthogonal property of the trial function to find 

spectral coefficients A . 

n 
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If the trial function Z^’s satisfy the following orthogonal i1 
relations 

b 

Xc3(x)Z(x)ZCx)dx = 6 (2.5.2) 

a 

where u(x) is a weight function and 6 is the Dirao-Delt 
function. Then by multiplying (2.5.1) by w(x)Z (x) (where 0<ni<N-l 
and integrating between the limits a to b we get 

b N-i b 

/ w(x) A(x) Z (x) dx = E S w(x) A Z (x) Z (x) dx 

m “ r» n Tw 

a a 

N-l 

= E A (C 6 ) 

ri r» nm 

n=o 

As m is between 0 to N-*l this sum will have only one non-zero ten 
when n=m. Therefore the sum will be equal to A C and 

r> n 

1 

K “ '■C?r~ ^ Z.^X) dx (2.5,3) 

a 

From equation (2.5.3) it can be seen that given a function A(x) anc 
a set of orthogonal trial function Z ’s. ¥e can find thf 
coefficients A^ for the approximate spectral expansion of A(x; b> 
integration. However, in the usual applications of spectral methods 
these integration will have to be done a very large number of 
times; hence to do them analytically is out of the question. 
Generally, therefore, a numerical quadrature rule is used for the 
integration in equation (2.5.3). 
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(2.5.b) Gaussian Quadrature 


The gaussian quadrature rules are the roost accurate in 
evaluating integrals by numerical integration .The accuracy of 
these method is a great jump over that of the regular interval 
quadratures because with N points we get 2N-1 polynomial accuracy. 
The idea of the gaussian quadrature is to get a good appro* i mat j on 
of an integral by weighted sum of the function values at certain 
specially chosen points i.e. 

b K-l 

f w(x) A(x) dx a: r w. ACx ) (2.5.4) 

, X x. 

a. V so 

where wCx) is a weight function of the orthogonal set Z (x), the 

TS 

A(x.) the function values at grid points x,’s and w *s are the 

». t. V 

weights assigned to those function values. 

The collocation points (also called the grid points) x ’s and 
weights w^ that give the most accurate evaluation of the integrals 
depend on the weight function co(x) and the limits a , b of the 
integral. Depending on the weight function used in the 

integral, we get various different quadratures like Gauss-Legendre , 
gauss- Chebyshev etc. and different corresponding collocation 
points . 

For an N point Gaussian quadrature the collocation points are 
chosen as the zeros of the order trial function Z^(x). It can 
be proved that these points get the highest order accuracy for a 
general integrand. In case of spherical harmonics(a trial function 
consisting of product of Legendre and Fourier terms), we have 
evenly spaced grids in the longitudinal direct ion(zeros of Fourier 
component) and unevenly spaced in meridional direction (zeros of 
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associated Legendre polynomial). 


The weight functions of the quadrature can be obtained by 
demanding that (2.5.4) be exact if A(x) is any of the orthogonal 
set Z <x),Z (x) ,Z (x). This means 

O l N-l 


o o o 


O A O 



+ w Z (x )= 

N-l O N-S. 

a 

w(x) 

Z dx 

o 


+ w„ ^ Z (x )= 

N-l t N-4 

Cl 

w(x) 

1 • • • • 1 

Z dx 

1 


wZ^ (x)+w Z (x)4 

O H- 1 O A N-1 1 


Z dx 

M- 1 M- i W- A <aL M- i 


Writing the above equations in matrix form ,we have 


[z] |w|=|cj (2.5.5) 

As x_, X ,x ^ are known (i.e. roots of the order 

O ft N— ft 

trial function ) the matrix IZl and column vector {C3 are also 
known. Thus by solving the matrix equation (2,5.5) ,we can easily 
determine weight functions of the quadrature i.e. {w}. 

Once the collocation points and weights are found for the 
quadrature to evaluate the integral (2.5.3), we can develop a 
relationship between function values at grid points and the 
spectral coefficients. 

Given that the values at the x. collocation point is A we see 

j J 

that the application of the Gaussian quadrature to (2.5.3) will 
yield 

J N - 1 

A = E w A. Z (X ) (2.5.6) 

n Cn . J J J 

J =o 
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which is called a "forward transform". Also given the values of the 
coefficients A we can evaluate the values at the grid points A by 

n J 

using (2.5.1) to get 

N-4 

A = E A Z (x.) (2.5.7) 

J n n t 

t\=o 


(2.5.c) Discrete Fourier Transform 

For a periodic function A(x) { 0 :£ x S L) there exists a 
finite transform called the discrete Fourier transform which relate 
the spectral coefficients Aj^ of the expansion 

N-i . , 2 , 

A^^^(x) = E A e^^ L ^ (2.5.8) 

k=o 

with the values of A(x) at the collocation points x^. These 
collocation points are defined as 

X. = (j L)/ N j=0,l ,N>-1 (2.5.9) 

The collocation points x^ are evenly spaced . 

If A^ is the set of functions values at the collocation points 

X. then 
J 


N- 1 


A. 


1 ^ „ -C2tTikj/N) 

TT~ E A. e 

^ j=o •' for k=0,l,2. 


(2.5.10) 


,N-1 


and the inverse transform is 


Aj = E A^ ^(2TTijk/N) 


(2.5.11) 
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Generally these summations in equations (2.5.10) and (2.5.11) 
are never done by a straight forward evaluation. The Fast Fourier 
Transform( an ingenious algorithm) is used for its evaluation. A 
straight forward evaluation takes O(N^) operations while a FFT 
takes only 0[N In^(N)] which is much smaller (for larger values of 
N). 

The FFT is a particular way of factorizing and rearranging the 
terms in the sums of the discrete Fourier transform. 

(2.6) .Aliasing 

In the spectral approximation generally the lowest-order 
coefficients are the largest and therefore most important. 

Therefore the larger the value of k less will be the significance 

of the corresponding coefficients. Given that two functions f^(x), 

g (x) are represented by their values at the grid or collocation 
points, how should their product h^(x)=f^(x) be evaluated ? 

The obvious answer is that the value of It^Cx) at the collocation 
points be obtained by multiplying the corresponding values of f„Cx) 
and g (x) at the same collocation points, i.e. 

)W 

h = f. g. , j= 1 N (2.6.1) 

J J -J 

However there is a subtle reason why this approach is not 
altogether right. The approximate functions and 

(N-l)th order polynomials and if their values at N collocation 
points is known the finite transform will exactly evaluate their 
corresponding coefficients. However if a polynomial of order 
greater than N-1 is represented on only N collocation points, the 
finite transform will no longer be able to evaluate the 
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coefficients exactly. All the coefficients will be have some error. 


If f^(x) and g^(x) are (N-1 )*■*"' order polynomials then 

obviously their product h^(x) will be of order 2N-2 . If h^(x) is 

represented by only N points there will be some error in the 
evaluation of the coefficients of b^(x) . This error is called 

alicteing error and effects both the low order important 
coefficients and the less important higher order coefficients. 

However, if f Cx) and g (x) are polynomials of order CN-l)/2 

each we can see that b^(x) will be polynomial of order N-1 and its 

coefficients can be evaluated exactly. This is done by finding 

coefficients from the collocation point values (f., g) and 

J J 

putting the upper half of the coefficients to zero to get the 
function values (f“, g“) at the collocation points. These function 

j j 

values will not be the original function values but will correspond 

to a polynomial of one-half the original order (because the upper 

half coefficients are put to zero). Therefore f* and g** are 

J J 

smoothened " values of the original values f . and g , and h. is 

J i J 

evaluated by 

h^ = f* g® , j= 1 ,2 N (2.6.2) 

As we know that the upper coefficients are the least 
significant and hence ignoring them will introduce the least error. 
In the original procedure of equation (2.6.1), the error affected 
all the coefficients of h^Cx) but in the second procedure (which is 
called truncation of modes) the lower half of the h (x) 
coefficients are exactly evaluated and the error is confined to the 
upper (less significant) coefficients. 
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(2.7) Spherical Harmonics 


For our purposes, spherical harmonics may be thought of as 
generalized Fourier Series applicable to the globe. The spherical 
harmonics have a number of useful properties that explain its wide 
spread use in the expansion of functions defined on a sphere. 
Namely, Spherical Harmonics (and its derivatives) are periodic in 
the longitudinal direction and are single valued and continuous at 
the poles. For infinitely smooth functions they also are spectral 
trial function i.e. the error in an expansion of a Spherical 
Harmonics drops off exponentially as the number of terms are 
increased . 


The spherical harmonic functions are functions of latit«de(<^) 
and 1 ongitude(X) . The order of the spherical harmonics is 
characterized by a pair of integers (n, m) . Spherical harmonics are 
so named because they are defined over the surface of a sphere and 
because solutions of Laplace’s equation were called Ha.x'momLc 
functions. Each Spherical Harmonic is the solution of the angular 
part of the Laplace equation. 


where 


Y (X,m) = P (P) 

r» 


i m X 
e « 


Mb sin4^ 


(2.7.2) 

(2.7.3) 


Spherical harmonics are therefore the product 
Legendre polynomial P"’(m) and the complex 

, * ^ 7 ’^ 

component) e'' a . 


of the associated 
exponential (Fourier 
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(2. 7. a) Associated Legendre Polynomials 


Because of their convenient properties for treating the 
latitude variation on a sphere, Associated Legendre Polynomials are 
often used in global atmospheric models. 


The associated 


n 



(2n +1) 
2 


legendre polynomials 
(n-m) I 

(n+m) 1 j 2^^ <n) ! 


may be represented by 


, m+ n 

~ m+n (2.7.3) 

d 


In this representation the functions are normalized as well as 
orthogonal i.e. 


4 . 

ptn pm d _ (2.7.4) 

r» k rs 

-1 

The function P"'(iu) has the following properties 

(1) There are n-m zeros of the function between the 
poles(/u=±l ) . 

( 2 ) Only this P“ are non-zero at the poles. 

(a) Those functions for which n+m is even are even 
(i.e. symmetric) and those for which n+m is odd are antisymmetric 
in (U. 


(2.7.b) Recursive relations for Spherical Harmonics 

For the the spherical harmonics Y*" = P"^ e'"™^ there are 

n n 

number of relations. In particular 
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0 


n < m 


(2.7 .13) 


_ 

ri 


Y“’^ = (-1)”' y"' * 

ri n 


C2.7. 14) 


n n 


(2*7. 15) 


(2.7. 16) 


Y s • 

n+l r»+i 


Y’^ _ y"* 
r\ ri-» n 


(2.7. 17) 


r 2n f 3 .. _ a. 

I 2n + 2 ) ^ ^ 


i/2 iX 


(2.7.18) 


The derivative of spherical harmonics along zonal direction is 
given by the relation 


if Yn 


I m Y 


(2.7.19) 


The derivative of Y^ along meridional direction is giveai by 

< 1- - 7 ^ — » < c. 

The laplacian of is a extremely simple relation 


7*Y"* 

n 


n (n-^l ) 

Z n 

a 


J2 

2 n 


(2,7-21) 


where = 


f _?_co.^+ -! 21— ] 

a* cos^ ^ d cos^i d 


(2.7.22) 
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It can be seen from (2.7.21) that Spherical Harmonics are eigen 
functions of the laplacian operator defined on a sphere. 

The orthogonality relationship for the spherical harmonics is 
given by 



Y” y" dp dX 

n j 


S 


( 2 , 7 . 23 ) 
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CHAPTER THREE 


NUMERICAL ALGORITHM 

This chapter presents the salient features of the numerical 
algorithm used in the present work. The topics covered include the 
spectral representation of U, V, 4*, and C fields, computation of 
the grid points, computation of the coefficients from the grid 
point values of the scalar fields and vice-versa, computation of 
the derivatives of the fields, the non-dimensional ization of the 
governing differential equations, and the spectral form of the 
governing differential equation. Back and forth transformation of 
4' and C fields, and computation of the Jacobian(non-l inear terms) 
are also brought out in this chapter. The later part of the 
chapter includes the time-stepping scheme, the flow chart, a 
description of the types of scalar spectra used in the analysis. 


(3.1) Spectral representation of scalar fields 

In this section, we present the spectral expansions of the 
scalar quantities used in our study. We use the so-called 
rhomboidal type of truncation. These variables are the steam 
functionCv) , the vorticity (C)» and the weighted x and y 
components of velocity (U and V respectively). The last two 
quantities are defined by 

U s u coa<P , V s V cos^ (3.1.1) 
where u, v are the components of the actual velocity, 4* is the 
latitude. The velocity components are so redefined, because (U, V) 
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are more convenient to use in the spectral barotropic equation. 
The spectral representation of the U, V, yi, and C fields are given 
below 


N- 1 N+ I m 


u = E 


n=|m 


ytn 


N- 1 K+ I m 

V = E Z 

m=-M+* n=lml 


y’« y’” 

n n 


N- 4 N+ Im I - 4 

V = E E 

m=-N+i n=lm| 


iT» -.rn 

yj y 


N-1 K+ Im J.-! 

C = E E , , C ^ n 

m=:->j+4 n=(ml 


(3.1.2) 


(3.1.3) 


(3. 1.4) 


(3.1.5) 


where Y*" = P"^(/j) e’’’"^ , the 5lpherioal Harmonic functions, ^ssin^, 

T% Y\ 

and U"*, v”* etc are the spectral coefficients of the corresponding 

n in 

physical quantities. In the above equations N is the rhoaboidal 
truncation number, and m and n are the 2-D wave numbers. The 
spectral coefficients are (complex) functions of time. The 
truncated series for U, V have one more term (in the n summation) 
than the series for w and C for the reasons that will become 
obvious later. 


(3.2) The Grid -point representation 

For a given truncation number N, some restrictions are 
imposed on the number of grid points used in the numerical model 
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in order to obtain alias-free computations. These restrictions are 

necessary for the non-aliased calculation of non-linear term 

Ce.g. .Jacobian term) using the transform method. For rhomboidai 

truncation N , the number of grid points in the longitudinal 

direction roust be greater or equal to 3N and in the latitudinal 

direction the number of grid points should be greater or equal to 
5N +1 

2 . These limits are built into the subroutine for the 

grid-spectral transforms. Therefore for N=16, we choose a grid 
point representation of 48x40Cthe first number corresponds to the 
longitudinal direction and the second in the latitudinal 
direction) . 


(3.3) Computation of the grid points 

For exact evaluation of a integral using Gaussian quadrature 
the grid points are chosen as zeros or the roots of an appropriate 
trial functions. In our problem C2-D on a sphere) the natural 
function to use are the spherical harmonics (eigen functions of 
poiseons equation on a sphere). The spherical tiarmonice are the 
product of an associated Legendre polynomial and a Fourier tor®. 
Therefore the grid points in the zonal or longitudinal direction 
are equally spaced (as for a Discrete Fourier transform), and in 
the meridional direction grid points are selected as zeros of the 
order Legendre polynomial (which are unevenly spaced) where £ 
is the numbpr of grid points in that direction. 
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(3.4) Computation of the coefficients from the grid points values 
and vice-versa 


The spectral representation of a function is 


ACV.P,) 


N- 1 


2 


N+ |m j - 1 

£ 

r,= jmj 


A P (u ) 

n n k 


V m X 

e J 


(3.4.1) 


defined on some suitable set of points K., and u . 

j k 

The spectral coefficients A*" are calculated from A(X,p), 
using the orthogonality property of the trial functions : 


, 271 i 

(3.4.2) 

Ti Z n n 

O - 1 

where the starred quantity is a complex conjugate. Here the limit 
of integration for X. varies from 0 to 27 t, and that for m from -1 
to 1. Generally the evaluation of the integral in (3.4.2) in 
carried out in two steps 


A’”(aj) = - J - ^ - S A(X,aj) dX (3.4.3) 

o 

and 

a"; = / A"'(p) P’"(m) d^J (3.4.4) 

n ^ 

- 1 

As we have function values at only a finite set of grid points, we 
can not directly use the aquations (3.4.3) and (3.4.4). Generally 
quadrature formulas (Gaussian, Trapezoidal , etc) are used to 
evaluate the above integrals. The equation (3.4.3) is calculated 
using a trapezoidal rule formula. Using this rule we can rewrite 
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equation (3.4.3) as 


= tV;V E * AC\ 5"” \ (3.4.5) 

j =o 

2 TT 

where = — g — j ,and j denotes the position of the grid points 

in the longitudinal direction and M is tne total aumber of grid 
points in the longitudinal direction and is equal to 3N. For any 
function which may be represented as a truncated Fourier series 
with wave number leas than M, this formula is exact. Generally, 
Fast Fourier Transforms iFFT) are used to evaluate (3.4.5). 


The Caxisss~L9gerudr& cj-xtadx'a.t\jix'» is used to evaluate equation 
(3.2.4). We can rewrite equation (3.2.4) using the Gaussian 
quadrature as 


A 


m 

n 


E w(^^) A'”(p^) p;;(p,^) 


(3.4.6) 


The sum is carried out over the K latitudes where are the 
roots of the ordinary Legendre polynomial of order K lie. 

=01 and are called Cauuseian lati txtdoe ; W(p^) are the Ccxxtssian 
weights. The Gaussian quadrature formula is exact for any 
polynomial of degree S 2K “1. The number of Gaussian latitudes are 
chosen as K=5M/2. It is clear that equation (3.4.6) is a 
multiplication involving a weight function matrixCwith elements 
V(/lJj^) 3, a coefficient matrixlwith elements A^'(M^) 1 , and an 
associated Legendre polynomial matrixlwith elements We 
have already explained the method of the computing weight function 
matrix in the chapter 2. We will later explain the comjiutat ion of 
the associated Legendre polynomial matrix or F matrix. 
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We can also compute the grid point values, if the spectral 
coefficients are known, by evaluating equation (3.4.1). Again 
this computation is performed in two steps. In the first step, we 
calculate using the relation given below 


a’-cp,) 


N+ I TH j - i 

£ 

|m| 



n 


f”’(R ) 
r» k 


(0,4.7) 


Equation (3.4.7) is basically a multiplication of a coefficient 
matrix and an associated Legendre polynomial matrix. The second 
step uses a transform method or FFT to find the grid point value 
of the function A(X.,u ). The equation used is given below 

J ^ 


ACX^.p^) 


W- 1 


£ 

Tn=-N-»- 1 


’(p,) 


(3.4.8) 


The transformation from the grid or the physical space to the 
spectral space is called the For-warci transform whereas that from 
the spectral space to the grid space is called or inverse 
transform. Generally, we call equation (3.4.5)-(3.4.6) as forward 
transform pair and equation (3 . 4 . 7)-(3 . 4 . 8 ) as backward transform 
pair . 


(3.5) Computation of the meridional and the zonal derivatives 

The reasons spectral methods have gained increasing 
popularity in applications requiring numerical solution of PDE’s 
is that they allow extremely accurate evaluation of derivatives. 
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Calculation of tho sonal d&rivatives 


Consider once again the spectral expansion 

ACX.p) = 5;: j: a"* 

ri ft 

m r^ 


(3.5.1) 


By differentiating (3.5.1) with respect to X., we get 


dA(X , <j) 


= E E 

m n 

= E E o’” 

“ “ ri TV 

m n 


( 3 . 5 . 2 :) 


Where d”= i m a” . From <3.5.1, 3.5.2) it is cl^ear that the 

spectral coefficients of the derivative (D”’) in the zonal 

r» 

direction are related to the original coefficients CA^') by a very 
simple relation. Once is found out the transform pair 

method [as in C3.4.7)-(3.4.8)] are used to calculate the values of 
the zonal derivatives nt grid points. 


Computation of th» m&ridional derivative 


By differentiating (3.5.1) with respect to p , we get 


g A ( A , ju ; ^ ^ ip n 

w ri 


(3.5.3) 


The derivative of an associated Legendre polynomial 
calculated using the recursive relationships. Once 

d P^Ai ^ 

derivative matrix, with element I)P^'(^<j^)= — g— , 


in (3.5.3) is 
we know this 

the transform 
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pair ( 3 . 4 . 7 )- ( 3 . 4 . 8 ) is used to calculate the grid point values of 
the meridional derivative of a given function. The computation of 
this derivative matrix or DP matrix of associated Legendre 
polynomial is explained in the next section. 


(3.6) Computation of Associated Legendre Polynomial and its 
Derivative 


The associated Legendre polynomial are a family of orthogonal 
polynomials in the interval (-1, 1) and have a weight 
function w = 1 . 


The associated legendre polynomials of the first kind 
normalized to unity are represented by 




(2ii +1) 
”12 


(n>m) I 

* J 2’^ ni dAi’"'"’" 




n 


(3.6.1) 


Associated Legendre polynomial (normalized one) satisfy the 
orthogonal relationship given below 

j. A p’w p« ^ (a. 6. 2 ) 

r> k nk 

The P”* rapidly become cumbersome to represent - the first few 

ri 

functions are 

p° = yu2 

0 

P° = ■/~372 iJ 

1 
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etc . 


p° = -rTrT (1/2) (3 (/- 1) 

It is clear from the above expressions that it will be btitter if 
we use some other means to generate the associated Legendre 
polynomials in a numerical model. Generally the P*” or P matrix 

T\ 

are all obtained using the following relationships 


P^ = 0 5 n < m 

P~” = (-!)•* p’" 

n f\ 

and the recursion relations 


(3.6.3) 

(3.6.4) 


(3.6.5) 

(3.6.6) 

] 

ii/z 


The elements of the P matrix used In the forward-backward 
transform are defined as P(m, n, ju ) s P^(k). The way this 
matrix is generated is as follows: For any given 

y i/a , then P° is computed using (3.6.5), and the higher n values 
for n=0 can be computed using (3.6.6) and the previously computed 
values of P. Then (3.6.7) can be used to generate and the 
process repeated using (3.6.5) and (3.6.6}. Similarly the P values 

at the m=2, 3, 4 levels can be generated. For the next the 

entire process will have to be repeated. 


(3.6.7 ) 

C3.6.8) 


VTITT-^ ^ P 


n+1 n+i ^ r\ n n***! 


rs+i 


where 


f 2 n_+ 

3 

t 2 n + 

T 

r- 2 2 


j n - m 


L 4 n*- 

1 
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Once we know the all p”(jUj^)’s, the derivatives of associated 
Legendre polynomials or DP matrix are calculated using the 
relation 


.. 7. dPn 

^ > -377- 


= - n £ p »(n+l) c'' p 

fi+i n+i n—i 


(3.6.9) 


dPn 


where -g - " at p is the element of the derivative matrix £>P(ii. n, 




In this way we can generate F and DP matrix .These matrix are 
used in backvxxrd. and forward transformation and in computing the 
meridional derivatives of the dependent variables. 


(3.7) Non-dimensional izat ion of the governing equation 


The problem under consideration in the present work is that 
of the non-divergent barotropic vorticity equation. The 
corresponding governing equation in spherical form is given below 


ar 

where 


1 

a cos*^ 



+ 


COSflb 


a 

a4> 




(3.7.1) 


COS0 ay 

a dtp 


(3.7.2) 


V = 


1 ay 

a a\ 


(3.7.3) 
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1 


d V 


(3.7.4) 


c * - ^ u 

0 ^ cos*^ 0 \ 

Here a is the mean radius of earth and O is the rotation rate of 
the earth. The primary directions are; longitude X increasing to 
the east, latitude ^ increasing to the north. 


For convenience equation (3.7.1) is changed to 
non-dimensional form by using radius of the earth as length scale 
and (2 n) as the time scale, as in Washington and Parkinson [ 18 ] . 
The quantities with * as superscript are non-dimensional 
quantities . 


* u * 

“To-a- ’ ^ = -T-OT 


* 


= W 


TIT 


> 


C 


* 


c 

nr 


Substitution in (3.7.1) ^ives 




at 


aoB^ip 


FjL 


(U 


* 

C ) 


cos^ 


0K 


¥e can write U and V in terms of the 
non-dimensional stream function : 


d 

a4> 


* « 1 

(VC )J- 


(3.7.5) 


derivatives of 
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As 


U = - 


COS<^ &H) 


dtp 


therefore , 

U* 2 O a = - 
which becomes 


-££5t- 2 n a’ 


U = - cos^ 




Otp 


and finally, 


U* e 


cos tp 


atp 




using 


a _ ajj a _ 

ap ap ap 


cosp 


whi le 


V = 


Sk 


becomes 


V* = 


a^> 

~Wk 


lit '*it 

substituting U , V .and 


4t-= coBp-^ in (3.1.5), we get 

ap 


»C 


at 


[- 


£1 ( ^ C*) + 


_£ ( ^ 


?*) ] - 


ay/ 

ax 


ax a/u dju ax 

Dropping ♦ for convenience, we can write above equation as 

aw 


at 


= -[ 


£_(_£iiL. c ) + 
ax a^j 


_a ( _a^ 

a^ ax 


C ) 


ax 


ac - _ 


(vj '> + 

a 

1 

i**-\ 

j 

1 - ^x 

at 

1 ax 


ap 

m 



On simp 


litication, the term inside the braokat can be rewritten as 
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f 


-^ = [ V). - ^ (3.7.6) 

Th6 term insiclB the bracket is called the Jacobian which is 
basically consists of the non-linear terms. 


(3.8) Governing equation in the Spectral form 

We now present the spectral form of the non-dimensional 
non-divergent barotropic vorticity equation used in the present 
study . 

We have already derived the non-dimensional form of the 
governing equation (3.7.6), and discussed the spectral 
representations of the scalar fields. We write the non-dimensional 
barotropic vorticity equation in the spectral form as 

ac’” 

—= F" - i m V'"’ , 05 « 5 N-1 (3.8.1) 

ja . Fi n 

^ ^ m<n<m+N-l 

where c"* and 1WN8 the complex coefficients of vorticity and 

ri r» 

stream function, respectively. The coefficients are the 

spectral coefficients of the Jacobian t * ^ ) 3 • Note 

that equation (3.8.1) need not be evaluated for negative m. This 
IS because we are dealing with real fields and this property can 
be used in finding, say, C~"' from by using 

r""’ = (- 1 )”' r"”* (3.8.2), 

n n 




(3.9) Computation of ip coefficients from C coefficients anc 
vice-versa 


The stream function v' is related to vortioity C by 

re 1 at ion 

C = V* V (3.9.1 

where C. V'l and the laplacian operator are al 
non-dimensional ized. Now 

7* = -n (n+1) Y"' (3.9.2 

ri n 

Considering (3.9.1) and (3.9.2) along with the spectra 
representations of ip and C I (3. 1.4) and (3.1.5)] we see that th 
spectral coefficients of ip and C are related by 

C’” = -n (n+1) (3.9.2) 

So given either v'"' or c’^ we can compute the other by multiply in, 

® n ri 

or dividing by -n (n+1), respectively. 


(3.10) Computation of the Jacobian 

The Jacobian comprises of the product of longitudinal am 
meridional derivatives of C and ip ) : 

JCC.V-) = Cv<^ c^) • <3.10.1) 
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The pseudo-spectral technique is used to evaluate this term. That 
IS, given the spectral coefficients of w and C. we first compute 
the grid point values of the longitudinal and meridional 
derivatives of these functions by using the techniques presented 
in section (3.5). The grid point values of the Jacobian are then 
evaluated from (3.10.1). The grid point values are then 
transformed to obtain the spectral coefficients of the Jacobian 
i.e. F^. This method is called pseudo-spectral because the 
multiplications are not done in spectral space but rather in 
physical (grid) space. 

(3.11) Flow-chart 

A flow chart given in Fig. (A) summarizes the algorithm used in 
the present work. We used pseudo-spectral technique in our problem 
(i.e., all multiplications have been performed on physical space 
or at grid points and all derivative calculations have been 
performed on spectral space). The time stepping scheme used is 
finite difference leap-frog scheme with a forward-difference 
start . 


(3.12) The spectrum of a scalar field 

One very important analytic tool is the study of unsteady 
velocity fields concerns the distribution of energy at the various 
length scales in the problem. The usual method for obtaining this 
information is to construct spectra of energy versus the wave 
number. These spectra can be three, two or one dimensional. In 
isotropic fluids, generally, one dimensional spectra are preferred 
because of convenience in plotting. In spherical coordinates two 
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]p = ■t'\°Tn€-s4-ejb 


pLOKi -'Ch^R I 







types of one-^d imens iona 1 spectra are defined naroely the so-callei 
spectrum which depicts the energy distribution in the wavi 

number of the coefficients and the 'n" spectrum, which does th< 
same for the ’n’ wave number. 


THm cp»ctT'\mts: 

Let H be a scalar field. We can write global average of H as 
< H > = — /f” J- ‘ H dp dX (3.12.] 


As H = X] £ H’” y"* , the equation (3.12.1) becomes 

" r» r» 
w n 


L H ^ H Y diu dX 

-1 rir» 

m n 


— E E s’” X * y”' dfj dX 

4n““n o -- - 

Tti n 


- 1 r» 


These integrals are zero for all m, n except m=0 and n=0, 




< H > = 


(3.12.5 


■/2 


¥e can also derive the expression for ensemble average of 
square of any scalar field H. ¥e write 


< H > = 


1 

'4 n 


27T 1 

X X 

o - 1 


H H dpi dX 


27T X ^ ^ 

tV E E E E h:; Y” dj. dx 

r» w p q O “1 


Using the orthogonality relationship we can rewrite this as 
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< H > = 


E E E E h" H*"* <5 6 

n m p q ^ r' n 


E E 

ri TO 


H H 

n r. 


(3.12.3 


where * indicates a complex conjugate. 


The fluctuating part of the scalar quant i ty IS 


H* = H - < H > . 




Substituting (3.12.2) and(3.12.3) in the above expression, we get 


gTO gTO* 

.-. < H** > = EE - (H° / 2) 

m n 


„TO „TO* 

... <h«> = ee-!v^ 


TO n 


(m=n^O) 


(3.12.4, 


The variance (or ’energy’) of the field is thus decomposed 


into its ’energy’ at individual wave numbers. This allows the 


-TO 


calculation of a 2-D spectrum E . 

n 



T.TO 


f h’'’ 

1 ^ 

gTO* 

r% 

for 

IB > 0 

n ^ |m i 


E 

Tt 

S= H 

[ (H"’ 

H"'*)/ 2 

ri 

f or 

m = 0 

ri > 0 

Now 

as £ IS 

ri 

the 

same 

as E”’^' for 

r> 

real 

fluids 

so we can sum 

al 1 

non-negative 

m to 

get 
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N-* |tn|4N-l 

^ H > «: 2 E (m=n=='0 excluded) 

Wk=0 -nsE I m j 

Similarly we can define the one dimensional ’m’spectrum 
that ^ 

< H*® > = E FCm) 

mso 


(3.12.5) 

F(m) such 

(3.12.6) 


where F(m) 


{ 


N-1 

E E° 

“ ri 
rn=l 

I Wl I +*J-l 


E 

n= [ml 


m = 0 

m > 1 


(3.12.7) 


We can also define the "n" spectrum G(n) such that 

t#2 


or 


2 2 

< > = E G(n) 

ri=o 


(3.12.8) 


where 


G (n)= i 


H 


N-l 


, .n=0 
z 


,in 


+ T I H I ,1 < n < N-l (3.12.9) 

n 

m=i 


N-l 

E 

Tn=n-M+i 


H 


m 


,N-1 < n < 2(N-1) 


The equations (3.12.7) and (3.12.9) are used to evaluate the 
spectra of of U, V, and C fields. The importance of m and n 
spectra lies in the fact that they gives us a good idea of the 
relative distribution of energy in the wave number space. 
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CHAPTER 4- 


RESULTS AND DISCUSSIONS 

In this chapter, we present results that demonstrate the 
effect of truncation of IC on the error growth, a comparison of 
resolution error and the initial condition dependent error <ICD), 
and a study of the relative rates of increase of the low wave 
number and the high wave number errors. 

Initial condition, in the form of say 42x41 spectral 
coefficients (per level) of u, v, P etc, are commonly obtained 
from objective analysis. Most of the energy in these fields is in 
the low wave number coefficients. ¥e seek to study the effect of 
"truncation" of these spectral coefficients on the growth rate of 
error. Truncation is the procedure by which only the spectral 
coefficients corresponding to a Rhomboidal truncation number N are 
used and the others are discarded. For example if vre truncate the 
IC at N=16, we will use only the lowest order 256 of the available 
42x41 coefficients. 

The truncation of IC causes the loss of information contained 
in the high wave number coefficients and hence introduces an 
initial error. ¥e truncate the "exact” IC at the rhomboidal 
truncation number N =16, and 32, to study the growth of error for 
the above two cases. The computations are performed on 192x160 
gr idCcorresponding to N=64 calculation). The vortioity field is 
found to be more sensitive to the truncation of IC compared to the 
velocity fields as there energy is distributed over a wide wave 
number range compared to the velocity fields which have there most 
of the energy in the low wave numbers only. 
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The CPU time increases very rapidly with an incretase in the 
resolution. As we know that our initial condition is most likely 
to be inaccurate in spite of our best effort, therefore there is a 
no point in increasing resolution beyond the point where initial 
condition dependent error (ICD) becomes dominant. A comparison of 
the growth of the resolution error and the initial condition 
dependent error(ICD) gives an insight into this problem. As our 
study shows that if the IC contains 10% error (Standard 
deviation=0. 1) , then it is not advisable to perform the 
computations on the 192x160 (N=64) grid because the computations 
on the 96x80 (N=32) grid also gives adequate resolution in 
simulation upto a period of three-days. To study the resolution 
error effect, we consider three sets of grids 192x160 <Na64), 
96x80 (N=32), and 48x40 (N=16). To study the ICD error, we 
artificially introduce errors (5% ,and 10% rms error) on "exact” 
IC and then simulate it for 3-days on 192x160 grid. 

A comparison of the rms error growth due to the errors in the 
low wave number coefficients and that in the high wave number 
coefficients shows that in the first case rate of error growth is 
faster than that in the second case. The sensitivity of the 
solution to the distribution of errorCor random number sets) is 
also carried out in the present study. The result shows that, 
although quantitatively there are slight variations, qualitatively 
the conclusions remains same. The method of perturbing the high 
wave number as well as the low wave number coefficients is 
explained later in this chapter. 


(4.1) Generation of Standard Initial Condition 

In this section, we present the method for generating the 
initial conditions for our problem. Our initial data was obtained 
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from objectively analyzed conditions for Feb 19, 1990 and wa# 
obtained from the NCMRWF, New Delhi. It contains 42x41 spectral 
coefficients per level for the u, and v fields. ¥e used the 500 mb 
level data for generating our initial conditions. Therefore we 
calculate the initial vorticity coefficients from these initial u 
and V coefficients. We use only 256 u, v coefficients 
corresponding to rhomboidal truncation number 16; these makeup 
about 95% of the total available energy in 500 mb velocity field. 
The reason for neglecting the higher modes of u, and v is tliat 
they will introduce unnecessary noise while calculating C 
coefficients from the u and v coefficients. This is because the 
vorticity depends on the differentials of the velocity components 
and the process of differentiation magnifies the importance of the 
the high wave wave number coefficients. Therefore the 256 u, v 
coefficients are used to determine the grid point values Con a 
48x40 grid) of U and V, and then of the meridional and zonal 
derivatives of U and V respectively. The grid point values of the 
vorticity are then computed using equation (3.7.4). These grid 
point values are then transformed to obtain the 256 vorticity 
coefficients . 

The C spectrum obtained with only 256 spectral coefficients 
is unphysical because it drops to zero after truncation wave 
number N=16. Therefore we artificially generate energy in the 
higher modes in the following manner; the 256 coefficients of C 
are used as initial condition and we run the program for 5-days on 
192x160 (N=64)grid. This gives us a smooth spectrum of vorticity. 
Also, as five days is the order of the large time-scales of the 
atmosphere we can expect that the high wave number regions of the 
spectrum to have reached an equilibrium by this time. This ensures 
that the spectrum obtained is physically valid in the higher wave 
numbers. If we had used 4096 coefficients of u, v to generate 4096 
coefficients of C, we would have obtained noisy unphysical 
coefficients in the higher wave numbers. In this way, we generated 
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4096 coefficients of C» which we will refer as the "exact” initial 
condition or the standard initial condition in the coming 
sections . 

It may be recalled, that for alias free computation the 
number of grid points should be 3N and 5N/ 2, in the longitudinal 
and meridional direction, respectively. So an N=64 computation can 
be done on a grid of 192x160 where the first number gives the 
number of grid points in the longitudinal direction and the second 
number in the meridional direction. 

Fig (la) presents plots of m spectrum of our initial U, V, f 
fields with N=64. FigClb) shows the m-apectrum of U on a 
Log-Linear scale at t=0 and t=3 days. 

(4.2) Method of measuring errors 

The error in computation (or prediction) of individual 
coefficients may occur either due to the propagation of error 
(like lack of observational data, inaccurate data etc.) in the 
initial data or due to bad resolution in the computation. 

To study the effect of the initial condition dependent (ICD) 
error on the accuracy of the prediction, we perturb the standard 
initial condition by a known amount; while the effect of 
resolution on the loss of predictability is studied using three 
different sets of grids e.g. , 48 x 40, 96 x 80, and 182 x 
ieO(corresponding to N=16, 32, and 64, respectively). For both 
cases, the error in the prediction of individual coefficient is 
calculated by defining 

s ( - f"" ) (4.2. i) 

n TH T\ 
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where is a coefficient obtainsd from the perturbed initial data 
or on the coarser grid and is a coefficient corresponding to 
the exact Solution obtained on the highest resolution and with 
the "exact” initial condition. Having computed the error 
coefficients for all (m,n) we can compute the ’error' spectrum 
in the usual way described in the last chapter. 

Now the ffrror variance is defined as 

< ( 7?"* - f"' )* (4.2.2) 

n,T« 

where 7 ? and ^ are the physical space variables whose coefficients 
are and ; the brackets < > indicate a space average. The 
error spectrum is related to the error variance in the manner 
shown in equation ( 3 . 12 . 5 ) in the last chapter. Fig(ic) shows the 
’m’ spectrum of error in U at t=0 and t=5 days for N=32 
computations. The error shown is the resolution error compared 
with the 'exact' (N=64) computation. 

Another measure of error used in the present study is a RMS 
(root mean square) error of the coefficients. The formula for 
calculation of RMS error is given below 

r » / Z ( T)"*- /£< Jf"’ >* (4.2.3) 

where r is a RMS error, and summations are done over all m, n. 


(4.3) Sensitivity to truncation of initial condition 


In this section, we investigate sensitivity of 
simulations to truncation of the initial condition. Our so 
exact initial condition contain 4096 spectral coefficients 
(corresponding to N=64). We consider here tv^o cases 
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truncation number N = 16 and 32) to study the effect of truncation 
of initial condition on the evolution of the model atmosphere. The 
computation performed usxng all the coefficients of our standard 
initial condition is used as a basis for comparison. Note that, ail 
the above computation are perrormed on 192x160 grid. Thus, cne 
resolution in all the simulations will be tiie same; tJie onii’ 
difference being in the initial condition error introduced by the 
truncation . 

In the first case, we use only inner 256 coefficients 
(corresponding to N =16) of our exact initial condition, but keep 
all the other parameters same. The 192x160 simulation has 4096 
spectral coefficients. Thus zero values are assigned to all the 
coefficients other than inner 256 coefficients of our so called 
exact initial condition. The computations are done for 3-days 
using this truncated initial condition. This simulation is 
compared with the "exact" simulation on the same grid using the 
full 4096 coefficients of the "exact" initial condition. The 
method of computing error spectrum and RMS value is already 
discussed in section (4.2). 

Similarly we repeat the above steps for the second case 
(N=32). The initial condition for tliis case is obtained by 
assigning zero values to all the coefficients other tJian inner 
1024 coefficients (corresponding to N=32) of exact initial 
condition 

Fig (2) present plots between RMS value and time for 
both cases. In the Fig(2) following notations are used : 

- (U/256), and (V/256) shows the results for U, and V 
fields using only 256 spectral coefficients of "exact" IC. 

- (U/1024), and (V/1024) shows the results for U, and V 
fields using only 1024 spectral coefficients of "exact" IC. 
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The RMS error for the U field due to the truncation of IC is 
17.2% and 5.2% for the two cases N=16, and 32, and for the V field 
it is 15.8% and 3.7% respectively. The corresponding figures for 
the vorticity is much higher, being 48.3% and 25.2%, respectively. 
The low RMS error for the U, and V fields compared to C field can 
be understood from the fact that the U, and V fields have their 
most of energies in the low wave number range. 

It is also clear from the Fig(2) that the RMS error curves 
are almost flat or they increase only slowly with respect to time. 
These result suggests that while truncation of initial condition 
(IC) introduces error in the computation, this error grows quite 
slowly with respect to time if sufficiently high resolution is 
used . 


(4.4) Comparison of Resolution error vs ICD Error 

In this section, a comparison is made betwieen the growth of 
initial condition dependent (ICD) error and resolution error. The 
ICD error is studied by deliberately introducing error of known 
RMS in our so called exact initial condition. ¥e consider the 5% 
and 10% (RMS) error case for this. We consider three different 
sets of grid (i.e. 48x40, 96x80,192x160) corresponding to N=i6, 32 
and 64, respectively, to study the resolution effect on growth of 
RMS error . 

Initial condition dependent (ICD) error is studied by 
deliberately introducing error on standard initial conditions. 
This error is introduced by randomly perturbing the coefficients 
of the "exact" initial condition so as to cause an error of a 
specified standard deviation from the original coefficients i.e. 
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^r, ~ ^ ^ ^ (4.4.1) 

where p is a pseudo random number of uniform distribution with 
zero mean and specified SDCstandard deviation), is a new 
initial spectra) co*>fficient with a perturbation about standard 
initial data, and is a coefficient of the standard initial 
condition of vorticity. We consider here two cases corresponding 
to SD( standard-deviation)= 0.05 and 0>10. This corresponds to the 
initial condition containing 5% rms error, and 10% rms error, 
respectively. We perform calculations on the 192x160 grid. 

Next we consider computations of resolution error. A study of 
this demonstrate the effect of grid— size on the accuracy of the 
solution of a given problem. As the number of grid-points is 
increased the accuracy of the solution increases, but the 
computing time increases very rapidly. Therefore an optimum 
selection of the grid which gives good accuracy with reasonable 
amounts of computation is very necessary. In these simulations, we 
consider computation on 192x160 grid as the "exact” one in 
comparison with simulations on coarser grid. This simulations can 
be used to compare the results obtained from computation on 48x40 
and 96x80 grids. The initial condition for 48x40 grid requires 256 
spectral coefficients, 1024 coefficients are required for 
computations on 96x80 grid and 4096 coefficients are required for 
192x160 grid. These coefficients are selected from standard 
initial data by truncation. The RMS value is calculated using 
(4.2.3). It is clear that even in the beginning of the simulation 
the RMS error for the truncated IC will not be zero. For example, 
in the 256 coefficient IC problem n™ and will be exactly same 
for the inner 256 coefficients but outside this region will be 
zero while will not be so. Thus the error will not be zero for 
all wave numbers and the RMS error will also be non-zero. 
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It is clear that, when the grid is made coarser, the IC also 
become truncated because a coarser grid can support less number of 
spectral coefficients. A question that immediately springs to mind 
is this : How important are the sjiectral coefficients that are 
discarded by truncation ? There are two aspects to the answer to 
this question. Firstly, it is true that the energy contained in 
the various wave numbers decrease dramatically with increase in 
the wave number. Therefore unless the truncation is at very low 
wave number (m<16), we would expect that the discarded 
coefficients to be of less importance than the retained ones. 
Secondly, we know that the higher order initial coefficients are 
more difficult to measure, because of various factors, and are 
therefore prone to be erroneous. Therefore when the higher order 
coefficients are discarded, not all of the loss will be real 
information. Therefore the ’worst case’ for truncated IC will be 
when all the discarded coefficients are ’exact’ (which is highly 
unlikely for real IC) and the ’best case’ will be when they are 
all erroneous. 

In keeping with the above arguments, we consider two cases in 
the resolution studies : The 'Worst* case scenario when all the 
4096 coefficients of the IC are exact. Thus every discarded mode 
contributes real loss of information. In this case the N=16 and 32 
simulations are compared with the N=e4 computations with the 4096 
coefficients of IC. In the ’Best’ case scenario we regard only the 
first 400 coefficients (corresponding to N=20) to be exact and all 
the rest to be completely erroneous. This assumption, though 
idealized, is closer to the situation in real IC. In this case the 
N=16 and 32 computations are compared to the N=64 computation with 
a 400 coefficient IC. There will be a truncation of IC, it may be 
noted only in the N=16 case (because the N=32 computation has 1024 
coefficients while the IC has only 400). This case is called "best 
case" because even the coarse grids in our study will capture most 
of the IC information as non-capturabl e large wave number range 
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have no energy. 


The results of the ICD and the resolution error studies are 
shown in fig (3a)-(3c), where the RMS error in U, V and C 5 
respectively, are shown versus time. The following notations ai*e 
used : 

- >»=16/(¥) and N=:je/(B) shows the results for the 

computations on the 48x40 grtd for the so called "worst" case and 
the "best" case respectively. 

- N=32/(:V) and M=32/(B) shows the results for the 

computations on the 92x80 grid for the so called "worst” case and 
the "best" case respectively. 

- N=64/0.05C¥) and N=64/0. 05(B) shows the results for 

the computations on the 192x160 grid when the IC contains 5% rms 
error for the so called "worst" case and the "best" case 

respectively. 

- N=64/0.10(¥) and N=64/0. lOCB) shows the results for 

the computations on the 192x160 grid when the IC contains 10% rms 
error for the so called "worst" case and the "best" case 

respectively. 


From Fig(3a) and figC3b) it is clear that the resolution 
error in U and V for the N=16 computation begins with an error 
12-18% .This error is larger than we would typically expect to be 
in the initial conditions. Also the initial error increases quite 
rapidly, reaching to 25-30% in three days. Therefore. the N=J 6 
coBputatioM offer degraded prediotJone both fron, the point of 
view of initial oondition error (due to truncation) and aUo due 
to poor resolution. However, when we observe the EMS value.s tor 
the N=32 simulations we see that the error level are lower than 
that in the simulation with lOX initial error for atleast three 
days even in the worst case. Indeed, for the best case =oenario 
the resolution errors in a N=32 calculations are below the ICD 
error for simulation beginning with 5* IC error for a period of 
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at least three days. This suggests that g 
5-10%, which is what we would expect in 
data, N=32 provides about the optimum 
predictions of upto three days. 


iven initial IC error of 
real objectively analyzed 
resolution required for 


However, these conclusions hold only for the velocity field 
not for the vorticity as figureOc) shows. It can be seen that 
even N-32 gives unacceptably high initial errors of 25% (worst 
case) and 10 % (Best case). This of course. is because the 
vorticity is very dependent on the small scale high wave number 
structure of the flow field and consequently contain much more 
energy in the higher wave number than does velocity field. This 
make is particularly sensitive to truncation which removes energy 
at the high wave numbers. 


This suggests that given about 10 % error level in the initial 
conditions we can perform simulations on an N=32 grid (92x80) 
rather than wasting computer time to do the same simulation on a 
f iner grid. 


Fig(3d)-(3f) shows the sensitivity of the IC error to random 
number sets (a) and (b). It can also be observed from the figures 
that the curves for the different random number sets are 
qualitatively similar. Therefore, above stated conclusions are 
independent of random sets used. 


(4.5) Low wave number vs High wave number errors 

In this section we investigate the wave-number dependence of 
the error growth. That is, first we introduce error only in the 
high wave number coefficients of the initial condition and then 
Q 0 ly j the low wave number coefficients and compare the growth of 
error in the two cases. Here we define 'Low wave number' 
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coefficients as those that corresponds to a truncation number of 
N=16. This includes the first 256 coefficients. ’High wave number’ 
coefficients are those that fall outside this region. In both 
cases we do two simulations each corresponding to rms error of 
approximately 5% on velocity fields. 


To introduce error in the high wave number or on the low wave 
number coefficients, we use pseudo random numbers of uniform 
distribution with zero mean and specified standard deviation. We 
have to calculate the standard deviation for the high wave number 
error (a ) case as well as for the low wave number error case 
which wiU give the specified SD (standard deviation) of 0.05 for 
the total field. The expression for calculating new SD is given 


as 


a * E„ 

H H 

1 


= a 


= a 


and 


(4.5.1) 


where E and E^ are the energy contained in the high and low wave 
numbers, respectively. E is the total energy m the initial 
condition (includes high wave number as well as low wave number 
coefficients). ^ is the standard deviation of the whole field 
(0.05) which we want to specify, and and ct^ are the SD for th • 
error introduced in the high wave number coefficients and the 
wave number coefficients respectively. The error is introduced on 
the desired coefficients using equation (4.4.1). In this way 
generate our new set of initial conditions. 

Again the computations are performed on the 192x160 grid for 
the 3~days. The rms error is computed using (4.3.3). 

Fig(4a) and Fig(4b) present plots of RMS error in U and 

V, respectively, vs time for the low wave nuhber (Lo) end ^ S 

^ 4. of the variation 

i-Tn 1 To also get an laea oi 

wave number (Hi) cases, lo oto s 


68 




FIG.(4a) HIGH-LOW WAVE No. RMS ERRORS(U-VE-:i 



tafc 


§ 3 S 3 

> > > > 
till 

5 5 ! S 



aouuB swy 



caused by different random number sequence (used in introducing 
the initial error), two different runs were performed for each 
case with with two different random number sequences (a) and <b). 
It can be seen from the figures that the initial error is in each 
case is not exactly 0.05. This is because error is introduced in 
the vorticity field (as we are using the vorticity equation) and 
the error in the U, V fields cannot be precisely determined a 
priori. However the initial errors are close enough for the 
purpose of qualitative comparison. 

It can be seen from the Fig(4a) and Fig(4b) that for the low 
wave number case the RMS error for both U and V fields increases 
much more rapidly than in the high wave number case. Therefore it 
can be concluded that, for approximately the same amount of total 
initial RMS error, the predictions are much more sensitive to the 
small wave number (large-scale) coefficients than the large-wave 
number (small-scale) coefficients. 


It can also be observed from the figures that the curves for 
the two different sets of random number (a) and (b) are 
qualitatively similar. Therefore, the above stated conclusions are 
independent of the random error introduced. 


CONCLUSIONS AND SCOPE FOR FURTHER WORK ; 


The present work has shown that the truncation number N=32 is 
an optimum resolution, if our IC contains 5-10% RMS error, for 
3-day prediction of the velocity field. However, the vorticity 
field requires higher resolution. It is also found that the 
predictions are much more sensitive to errors in the lower order 
coefficients than in the higher order ones. The error introduced 
by truncation of IC is found to have slow growth rate if 
sufficiently high resolution is used. 
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The random number dependence test has also been carried out in 
our study. We observec that the results are qualitatively similar 
and are, thus, independent of the random number sets used to 
introduce error on IC. 

A very interesting extension to the present work would be the 
study of error spectra, the distribution and growth of error in 
the wave number space. During the course of the present work, this 
study has been investigated to a limited extent. It would also be 
interesting to do similar resolution v/s ICD error comparison 
using more realistic models. 


70 



BIBLIOGRAPHY 

1. Ttiompson, P.D., 1957: LfricertaiTity of zniital stccte cxs a footoT 
i,Ti tfxB pj'odtc iaht I i i y of 1 cxrg& scat& czimosph&X'i c: flow pat toms. 

Tellus, 9, 275-295. 

2. Loenz, E.N. , 1969: Th» pi'&d.ictajbilxty of a flaw ch pc:>$%'B€?ss»&SB 
many scales of motion, Tellus, 17, 321-333. 

3. Kraichnan, R.H. , 1970: InstatiiLi ty in fully develop&ci 

turbultfnc^ , Ph.ys. Fluids, 13, 569-575. 

4. Leith, C.E., 1971: AtmDsphBrz.c predic taJbi li ty and tv>o 

din^nsional turbulence, J. Atmos. Sci., 28. 145-161. 

5. Leith, C.E. , and R. H. Kraichnan, 1S72: Predictability of the 
txxrbulent flows, J. Atmos. Sci., 29, 1041-1058. 

6.. Basdevant , C. , B. Legras, R. Sadourny and M. Beland 1981: A 
study of the Barotropic model flows: Intermi t tarioy , xoaves 

and predictability, Atmos . Sci . , 38, 2305-2326. 

7. Holloway, G. , 1983: Effects of the planetary waue propagation 
and finite depth in the predictability of atmospheres. J. Atmos. 
Sci., 39, 314-327. 

8. Vallis, G.K. , 1983: On the predictability of guasz.~seosiroph^c 

flow: The effect of beta and baroclinicity , J. Atmos. Sci. , 40, 

10-27. 

9. Roads, J.O. , 1985: Temporal uariations in predz.c tabi li ty , J. 

Atmos. Sci., 42, 884-903. 

10. Strauss, D.M. , 1989: Baroclznic instability and waue-waue 

interactions in quasi— geos trophic error growth, J. Atmos. Sci., 


71 



42, 2380-2403. 


11. Charney, J.G. , R.G. Fleagle, V.E. Lalley, H. Riehl and D.Q. 

Wark , 1966: T?ve feasibility of a global observational and and 

analysis experiment, Bull. Amer. Meteor. Soc., 47, 200-220. 

12. Smagoriusky, J.,1969: Problems and promises of dei ermini si i c 
extended, range forecas t ing , Bull . Aoier. Soc . , 50, 286 — 311. 

13. Williamson, D.L. , and A. Kasahara. 1971: Adoption of 

meteorological varidhles forced by vpdatxng , J . AtmOS . Sci . . 28, 

1313-1324. 

14. Baede, P. , D. Dent and A. Hoi 1 ingworth, 1977: The effect of 
metric precision on some Tneteoro logical integration, 

15. Orszag, S.A,, 1970: Transform method for the calcvlation of 

vector-coupled sums: application to i/i© spectral form of the 

vorticity equation, J. Atmos. Sci., 27, 890—895. 

16. Haltiner, G.J., and R. J. Williams, 1980: Numerical Prediction 
and Dynamic Meteorology , John Wiley and Sons, New York. 

17 Holton, J.R. , 1979: Introduction to Dynamic Meteorology, 

Academic Press, New York, 

18. Washington, W.M., and C.L, Parkinson. 1986: An Introduction to 
Three-Dimensional Climate modeling. University Science Boohs, 
Oxford University Press. 


72 



